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In this paper we consider the anharmonic corrections to the anisotropic elastic rod model for 
DNA. Our model accounts for the difference between the bending energies of positive and negative 
rolls, which comes from the asymmetric structure of the DNA molecule. We will show that the 
model can explain the high flexibility of DNA at small length scales, as well as kink formation at 
high deformation limit. 



Characterizing the elastic behavior of DNA molecule 
is of crucial importance in understanding its biological 
functions. In recent years, single-molecule experiments 
such as DNA stretching and cyclization P, 0] have pro- 
vided us with valuable information about the elasticity 
of long DNA molecules. The results of these experiments 
can be described by the elastic rod model (also called 
wormlike-chain model) 0, 01 • In this model it is assumed 
that the elastic energy is a harmonic function of the de- 
formation 0, 01 . The elastic rod model is very successful 
in explaining the elastic behavior of the micron-size DNA 
molecules. 

Recently, modern experimental techniques have made 
it possible to study the elasticity of DNA at nanometer 
length scale [1, 0,11, 0- In these experiments it is ob- 
served that short DNA molecules are much more flexible 
than predicted by the elastic rod model. Several different 
models have been presented by now, that try to explain 
the origin of this discrepancy by considering the possibil- 
ity of local DNA melting [^, [l^, , or the occurrence 
of kinks in the DNA structure [l3| . Also Wiggins et al. 
have suggested an alternative form for the elastic energy 

Since the DNA is not a symmetric molecule, the en- 
ergy required to bend the DNA over its major groove 
is not equal to the energy required to bend it over its 
minor groove. The model which in introduced in this 
letter takes this difference into account. The effect of 
asymmetric structure of DNA on its elastic energy has 
been discussed previously by Marko and Siggia Q , where 
they showed that there must be a coupling term between 
bend and twist in the harmonic elastic energy. We will 
discuss that the asymmetric structure of DNA can also 
be introduced as a correction to the harmonic elastic en- 
ergy, which is of the third order. We shall show that our 
asymmetric elastic rod model can account for the high 
flexibility of short DNA molecules. 

In the elastic rod model DNA is represented by a con- 
tinuous inextensible rod. The curve which passes through 
the rod center determines the configuration of the DNA 
in three dimensional space. This curve is denoted by r, 
and is parameterized by the arc length parameter s (Fig- 
ure [T|) . In addition, a local coordinate system with an 
orthonormal basis {di,d2,d3} is attached to each point 



of the rod. As depicted in Figure [1] ^3(5) is tangent to 
the curve r at each point. c?3(s) = dr/ds, c?i(s) is per- 
pendicular to (^3(5) and points toward the major groove, 
and d2{s) is defined as ^2(5) = ^3(5) x di{s). These three 
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FIG. 1: Parameterization of the elastic rod. The local frame 
{di, d2, da} is attached to the rod. 

orthogonal vectors uniquely determine the three dimen- 
sional configuration of DNA. From classical mechanics 
we have 



CI X di 



1,2,3, 



(1) 



where the dot denotes the derivative with respect to s, 
and a is called the spatial angular velocity. The com- 
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ponents of f2 in the local coordinate system are denoted 
by f2i, fl2, and Ct^. The elastic energy of an inextensible 
DNA in the most general form can then be written as 



E : 



£[r!i(s),r!2(s),f^3(s)]ds, 



(2) 



where L is the total length of DNA. £{s) is the energy per 
unit length of DNA, i.e the energy density, at point s. For 
small deformations, the energy density can be written as 
a Taylor expansion about the lowest energy configuration 
. For a DNA with no intrinsic curvature and a constant 
intrinsic twist luq , the lowest energy configuration is given 
by do = [0,0,a;o]'^- Thus, at the lowest order, we arrive 
at a harmonic energy density in the form 
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(3) 



where fc^ is the Boltzmann constant, T is the tempera- 
ture, and w is defined hy tu = Cl — Qq. Qisa3x3 sym- 
metric matrix whose elements are the elastic constants 
of DNA 0, i, Considering a short segment of DNA 
with the length ds at the point s, this segment has a 
symmetry under 180° rotation about the local di axis at 
the point s. Thus the odd powers of fli must not appear 
in the expansion of energy density, and the matrix Q has 
only four independent non-zero elements: Qn, Q22, Qasj 
and Q23- Therefore, the harmonic energy density can be 
written as 



arm = -^kBT[Ai 17? + ilj + ^(^(3 - iJof 



-2Dn2{n3-ujo)]. 



(4) 



The first two terms in equation Q correspond to the 
bending of DNA over its grooves (roll), and over its back- 
bone (tilt), respectively. Ai and A2 are the correspond- 
ing bending constants. Since roll requires less energy 



121, ll8|, one expects that A2 < Ai. The 



than tilt [ 

third term indicates the energy needed for twisting the 
DNA about its central axis, with the twist constant C. 
Finally, the fourth term accounts for the coupling be- 
tween roll and twist [l^. Although the elastic constants 
of DNA may depend on the sequence [l^ , in this paper 
we neglect sequence dependence, and assume that they 
are constant all along the DNA. 

The existence of twist-roll coupling indicates that there 
is indeed a difference between bending over major groove 
(positive roll), and bending over minor groove (negative 
roll): For positive values of D, the DNA has a tendency 
to untwist when roll is positive, and to overtwist when 
roll is negative. 

To account for the effect of asymmetry on the bending 
energy of DNA, we need a term in the energy density 
which is an odd function of n2, and does not depend 
on r^i or r^a. There is no such term in the harmonic 
elastic energy, so we consider third-order terms in the 



expansion of energy density. The term proportional to 
ilf has the desired property. On the basis of some the- 
oretical analysis [20|, as well as experimental evidences 
21 1 and simulation studies [2^ , we assume that negative 
roll is more favorable than positive roll. Thus we write 
the third-order term in the form +1/31 where F 

is a real parameter. (It must be noted that the main 
conclusion of the paper remains valid if positive roll is 
easier than negative roll. To account for this case, one 
can write the third-order term in the form — 1/3! llf-) 
To keep the model as simple as possible, we neglect cou- 
plings in all orders, as well as higher-order corrections to 
the twist energy. So the only third-order term which en- 
ters in the model is 1/3! fij. Since the elastic energy 
must have a lower bound, we must keep the fourth-order 
correction to the roll energy, i. e. the term proportional 
to riji in the model. For consistency of the model, we 
also keep the corresponding fourth-order correction to 
the tilt energy. Since the anisotropy in bending energy 
is accounted for in the second order, to reduce the model 
free parameters, we write the fourth-order terms in the 
form 1/4! G'^(rif -I- rij), with G real and positive. Adding 
third-order and fourth-order terms to the harmonic en- 
ergy density, we obtain the asymmetric elastic rod model 
which is given by 
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This model accounts for the asymmetry between posi- 
tive and negative rolls, as well as the difference in the 
energies of roll and tilt. Since there is no coupling term 
in the model, roll, tilt, and twist can be regarded as in- 
dependent deformations, and the energy density can be 
decompose into three separate terms 



*.^asym 

where 



£1 ~ ksT 



rl 



1 



-A,ni + -G'nt 



£2 = knT 



;:A2nl+ iF'nl + ^G'nl 



1 

4! 



£3 = i^ksTCin^ - ujq) 



(6) 



(7) 



(8) 



(9) 



Depending on the values of ^2, F and G, £2 can take 
three different forms (see Figure [2]). For small values of 
F, £2 has only one minimum at f22 = and its curvature 
is always positive. For F > (2 j42 G'^)"'^/^ the curvature 
of £2 can change sign and there are two inflection points. 
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For given A2 and G there exists an upper bound Fc = 
(8/3 A2G^y/'^, beyond which £2 has two minima: one 
at = and the other one at a negative U2- In this case 
DNA has two stable configurations, and there is always 
a barrier between them. However, one expects that the 
barrier is not large compared to fc^T for a real DNA. 




FIG. 2: (Color online) £2 as a function of 0,2 for A2 — 
43.50 nm and different values of F and G. Black: F = G — 0. 
Blue: G — 3.20 nm and F = 7.20 nm, £2 has one minimum 
and its curvature is positive every where. Green: G = 3.20 nm 
and F = 7.80 nm, £2 has two inflection points and one min- 
imum. Red: G = 3.20 nm and F = 7.90 nm, £2 has two 
minima. Magenta: G — 3.20 nm and F — 7.94 nm. Note that 
a change of 0.04 nm in F results in a change of about 1 kBT/lo 
in £2 in the vicinity of the second minimum. 

To study the elastic behavior of DNA in the asym- 
metric elastic rod model, we calculate the distribution 
function P{9), the probability that the DNA bends into 
an angle 9. We use a Monte Carlo simulation to cal- 
culate P{9). In this simulation we discretize each chain 
into separate segments of length Iq = 0.34 nm, equal to 
the base-pair separation in DNA. The orientation of each 
segment is then determined by a vector 0, where |8| de- 
termines the rotation angle of the local coordinate sys- 
tem with respect to the laboratory coordinate system, 
and the direction of indicates the normal to the plane 
of rotation. The special angular velocity is related to Q 
as fl = Q/lo- In each Monte Carlo move, we randomly 
choose a segment along the chain, and for that segment 
we change the vector Q by A6. The direction of AO 
is random, and its magnitude is chosen randomly in the 
interval [0 , Qq] . Qq is chosen so that the accept ratio 
is about 0.5. We do not include the self avoiding in the 
simulation, since the probability of self crossing is small 
for the short simulated DNA molecules. 

Recently, Wiggins et al. have used atomic force mi- 
croscopy to measure distribution of the bending angle of 
short DNA molecules Q. Although the DNA molecules 



in Wiggins et al. experiment has the characteristic prop- 
erties of two dimensional polymers [1], to simulate the 
experiment we do not confine the DNA completely in a 
plane. The reason is that the minimum energy configu- 
ration of an anisotropic DNA is not plannar, although 
the deviation from a plannar configuration is negligi- 
ble [2^. It is known in the anisotropic harmonic elas- 
tic rod model, that the efi^ective bending constant of 
long DNA molecules in three dimensions is equal to the 
harmonic average of Ai and A2 [H, 17, ^off = 
[1/2 (1/^1 + 1/^2)]"^, while in the two dimensions the 
effective bending constant is given by A^^ = \J A\ A2 
[2^ . Since is always greater than ^eff, confining the 
DNA in a plane costs energy. For this reason, we allow 
the DNA to come out of the plane by 0.3 nm, which is 
seven times smaller than DNA diameter and lies in the 
range of atomic length scales. 

Following other studies [i, 25 1, we assume A\ ~ 2A2, 
C = 100 nm, and ujq = 1.8 nm^^. The values of A2, 
F and G are then determined by fitting the theoretical 
results to the experimental data of Wiggins et at, with 
the constrain that the persistence length of the DNA is 
54 nm . Figure [3] shows a good fit to the experimental 
data, which corresponds to A2 = 43.50 nm G = 3.20 nm 
F = 7.90 nm. The predictions of isotropic-harmonic elas- 
tic rod model are also shown in the Figure. We report the 
values of the fitting parameter with three significant dig- 
its. The reason is that £2 is very sensitive to the changes 
of the parameters, specially when it has two minima. In 
fact, a change in the order of 10~^ nm in these parameters 
may results in a change in £2 in the order of 1 ksT/lo (see 
Figure [2|), and therefore can significantly affect the elas- 
tic behavior of DNA. We must note here that the ratio 
A1/A2 is also relevant to the fitting procedure. However, 
one can obtain equally good fits for different values of 
A,/A2. 

The functional form of £2 for A2 = 43.50 nm, F ~ 
7.90 nm and G = 3.20 nm is shown in Figure [21 which 
has two minima. The second minimum occurs at ^2 = 
— 3.3nm~^ which corresponds to a —64° roll between 
adjacent base pairs. Thus, the existence of a second 
minimum can lead to the formation of kinks in the mi- 
nor groove direction in a tightly bent DNA. As can be 
seen in Figure [3l both in the experiment and our model, 
— log P(6') deviates from the harmonic model at large 
bending angles. Continuing the graphs in our model, 
they arrive to an approximately linear regime. This lin- 
ear behavior is a signature of kink formation. Both the 
slope of the line, and the crossover point are related to 
the values of F and G. 

The possibility of kink formation in the DNA structure 
has been considered previously by other authors. It was 
firstly mentioned by Crick and Klug, who proposed an 
atomistic structure for a kinked DNA [23|. They sug- 
gest that DNA can be kink most easily toward the minor 
groove. Nelson, Wiggins, and Phillips have presented a 
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FIG. 3: (Color online) — log P(6) as a function of 6, for L — 
5nm (red), L = 10 nm (green), and L = 30 nm (blue). Dots 
show the experimental data of Wiggins et al. . Curves show 
the theoretical results. Dashed curves: isotropic harmonic 
elastic rod model, with Ai — A2 — 54 nm, and F — G — 0. 
Solid curves: Ai = 87.00 nm, A2 = 43.50 nm F = 7.90 nm, 
and G — 3.20 nm. 



simple model for kinkable elastic rods [13[, in which the 
kinks are completely flexible, and can be formed in any 
direction with equal probability. Their model can explain 
the high cyclization probability of short DNA molecules 

f. The linear behavior is also observed in their model 
but the slope of the line is always zero. In a re- 
cent experiment, Du et al. proved the existence of kinks 
in DNA minicircles of 64-65 bp [2^. Molecular dynam- 
ics simulations on a 94 bp minicircle [13] also show that 
kinks are formed, with the same structure predicted by 
Crick and Klug. Similar kinks in the direction of minor 
groove have been observed in the structure of nucleoso- 
mal DNA [U. 

Kinks are also observed in the cr ysta l structures of non- 
histone protein-DNA complexes [l9l . [27j . In these com- 
plexes, DNA has a clear tendency to kink in the ma- 
jor groove direction. Du et al. have found the distri- 
bution function P{9) for a base-pair step in these com- 
plexes [28|. Although it is contradictory to our primary 
assumption, that kinks are formed toward the major 
groove, the asymmetric elastic rod model can be fitted 
to the Du et al. data by writing the third-order term 
in the form — l/SIF^Slj, and choosing F = 8.86 nm and 
G = 3.83 nm. We found that the model with these val- 
ues of F and G can not explain the experimental data of 
Wiggins et ai. This shows that the statistical property 
of DNA in the protein-DNA complexes certainly differs 
from the free DNA. This difference is probably due to 
the interaction of proteins with DNA, which can alter 
the DNA conformation dramatically, and leads to differ- 
ent effective elastic constants. 



In this paper, we presented a generalization of the 
anisotropic elastic rod model, assuming that the ener- 
gies of positive and negative rolls are different as a result 
of the asymmetric structure of DNA. We showed that 
this model can explain the elastic behavior of short DNA 
molecules. We also showed that this model allows the for- 
mation of kinks in the DNA structure when the DNA is 
tightly bent. The kinks always form in one of the groove 
directions, as suggested by other studies. 
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